#coding=utf-8
# ////相关分析
import numpy as np
import arcpy
# 降水栅格数组
rainList=[]
# ndvi栅格数组
ndviList=[]
for x in range(15):
    rainrasterpath='C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\temRaster\\'+str(2005+x)+'data.tif'
    rainimage=arcpy.RasterToNumPyArray(rainrasterpath)
    rainList.append(rainimage)
    ndviimagepath='C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\15yearndvi\\'+str(2005+x)+'ndvi.tif'
    ndviimage=arcpy.RasterToNumPyArray(ndviimagepath)
    ndviList.append(ndviimage)
# ////相关分析
# 求影像1和影像2间的相关系数 first 为降水栅格或者温度栅格，second为ndvi栅格
def corrToraster(imageList1,imageList2):
    first=imageList1[0]
    row=first.shape[0]
    column=first.shape[1]
    # 对每一个像素进行操作
    # 初始化一个栅格图
    initRaster=np.zeros( (row,column) )
    for rowindex in range(row):
        for columnindex in range(column):
            # 如果检测到nodata区域则不计算这个位置 
            if(first[rowindex][columnindex]<-200):
                    continue
            fisrtarr=[]
            for image in imageList1:
                fisrtarr.append(image[rowindex,columnindex])
            secondarr=[]
            for image in imageList2:
                secondarr.append(image[rowindex][columnindex])
            x=np.vstack((fisrtarr,secondarr))
            # //得到同一个像元上温度或者降水与ndvi的相关系数阵
            re=np.corrcoef(x)
            # 将计算得到相关系数赋给新栅格
            initRaster[rowindex][columnindex]=re[1][0]
    myRaster = arcpy.NumPyArrayToRaster(initRaster)
    myRaster.save("C:/Users/xulzu/Desktop/paper/origin_data/iamge/temNdviscale.tif")
corrToraster(rainList,ndviList)



# % nullbloack=[0,1,2,3,4,5,9,10,11,12,13,14,20,21,22,23,30,31,40,68,69,77,78,79,87,88,89,95,96,97,98,99];
# 最后一步运算

import numpy as np
import arcpy
# 计算NDVI和时间的相关系数
def polyfit(rainList,block):
    # 获取影像的行列数目
    first=rainList[0]
    row=first.shape[0]
    column=first.shape[1]
    # 初始化一个二维数组
    initRaster=np.zeros((row,column))
    # initIntercepterRaster=np.zeros((row,column))
    # 对每一个栅格单元进行迭代以计算相关系数
    for rowindex in range(row):
        for columnindex in range(column):
           # 如果检测到nodata区域则不计算这个位置 
            if(first[rowindex][columnindex]<-200):
                    continue
            fisrtarr=[]
            # 读取每一年中栅格单元ij上的值形成一个数据列表
            for image in rainList:
                fisrtarr.append(image[rowindex,columnindex])
            # 整理出时间类表
            secondarr=[2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019]
            # for image in ndviList:
            #     secondarr.append(image[rowindex][columnindex])
            # cor=np.polyfit(fisrtarr,secondarr,1)
            # initscaleRaster[rowindex][columnindex]=cor[0]
            # initIntercepterRaster[rowindex][columnindex]=cor[1]
            # 下面开始计算相关系数
            x=np.vstack((fisrtarr,secondarr))
            # //得到同一个像元上温度或者降水与ndvi的相关系数阵
            re=np.corrcoef(x)
            # 将计算得到相关系数赋给新栅格
            initRaster[rowindex][columnindex]=re[1][0]
    myScaleRaster = arcpy.NumPyArrayToRaster(initRaster)
    myScaleRaster.save("C:\\Users\\xulzu\\Desktop\\gismiddata\\15ndviadd\\ndviAddRaster"+str(block)+".tif")
    # myIntercepter =arcpy.NumPyArrayToRaster(initIntercepterRaster)
    # myIntercepter.save("C:\\Users\\xulzu\\Desktop\\gismiddata\\interceptRaster\\interceptRaster.tif")

alliamgeList=[0,1,2,3,4,5,9,10,11,12,13,14,20,21,22,23,30,31,40,68,69,77,78,79,87,88,89,95,96,97,98,99]
start=90
循环对每一份切割后的影像计算相关系数
for index in range(start,start+10):
    if index in alliamgeList:
        continue
    # 降水栅格数组
    rainList=[]
    # ndvi栅格数组
    ndviList=[]
    for x in range(15):
        rainrasterpath='C:\\Users\\xulzu\\Desktop\\gismiddata\\clipdata\\15yearclip\\'+str(2005+x)+'data'+str(index)+'.tif'
        rainimage=arcpy.RasterToNumPyArray(rainrasterpath)
        rainList.append(rainimage)
        # ndviimagepath='C:\\Users\\xulzu\\Desktop\\gismiddata\\clipdata\\ndviRaster\\'+str(2005+x)+str(index)+'ndvi.tif'
        # ndviimage=arcpy.RasterToNumPyArray(ndviimagepath)
        # ndviList.append(ndviimage)
    polyfit(rainList,index)


initscaleRaster=[]
# 影像拼接
alliamgeList=[0,1,2,3,4,5,9,10,11,12,13,14,20,21,22,23,30,31,40,68,69,77,78,79,87,88,89,95,96,97,98,99] 
zeroarr=np.zeros((543,683)) 
columnarr=np.zeros((0,6830))
for row in range(10):
    rowarr=np.zeros((543,0))
    for column in range(10):
        index=row*10+column
        if index in alliamgeList:
            rowarr=np.hstack((rowarr,zeroarr))
            continue
        iamgepath='C:\\Users\\xulzu\\Desktop\\gismiddata\\15ndviadd\\ndviAddRaster'+str(index)+'.tif'
        readimage=arcpy.RasterToNumPyArray(iamgepath)
        rowarr=np.hstack((rowarr,readimage))
    columnarr=np.vstack((rowarr,columnarr))
myScaleRaster = arcpy.NumPyArrayToRaster(columnarr)
myScaleRaster.save("C:\\Users\\xulzu\\Desktop\\gismiddata\\ndviAddRaster.tif")

for index in range(15):
    ndvipath='C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\15yearndvi\\'+str(2005+index)+'ndvi.tif'
    rainpath='C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\rainRaster\\'+str(2005+index)+'data.tif'
    result=Raster(ndvipath)-Raster(rainpath)*Raster('C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\iamge\\rainCorRaster.tif')-Raster('C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\iamge\\interceptCorRaster.tif')
    result.save('C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\15yearsResidual\\'+str(2005+index)+'data.tif')


# 人类活动一元线性回归
import numpy as np
import arcpy
# ////相关分析
# 求影像1和影像2间的相关系数 first 为降水栅格或者温度栅格，second为ndvi栅格
def corrToraster(imageList1):
    first=imageList1[0]
    row=first.shape[0]
    column=first.shape[1]
    # 对每一个像素进行操作
    # 初始化一个栅格图
    initRaster=np.zeros((row,column))
    for rowindex in range(row):
        for columnindex in range(column):
            # 如果检测到nodata区域则不计算这个位置 
            if(first[rowindex][columnindex]<-200):
                    continue
            fisrtarr=[]
            for image in imageList1:
                fisrtarr.append(image[rowindex,columnindex])
            secondarr=[2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018,2019]
            # cor=np.polyfit(fisrtarr,secondarr,1)
            # initRaster[rowindex][columnindex]=cor[0]
            x=np.vstack((fisrtarr,secondarr))
            # //得到同一个像元上温度或者降水与ndvi的相关系数阵
            re=np.corrcoef(x)
            # 将计算得到相关系数赋给新栅格
            initRaster[rowindex][columnindex]=re[1][0]
    myRaster = arcpy.NumPyArrayToRaster(initRaster)
    myRaster.save("C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\iamge\\residualScale.tif")
# ndvi栅格数组
ndviList=[]
# rastersize=[]
for x in range(15):
    ndviimagepath='C:\\Users\\xulzu\\Desktop\\paper\\origin_data\\15yearsResidual\\'+str(2005+x)+'data.tif'
    ndviimage=arcpy.RasterToNumPyArray(ndviimagepath)
    # rastersize.append(ndviimage.shape)
    ndviList.append(ndviimage)
corrToraster(ndviList)

